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Abstract. Self-similar processes are useful in modeling diverse phenomena that exhibit 
scaling properties. Operator scaling allows a different scale factor in each coordinate. 
This paper develops practical methods for modeling and simulating stochastic processes 
with operator scaling. A simulation method for operator stable Levy processes is devel- 
oped, based on a series representation, along with a Gaussian approximation of the small 
jumps. Several examples are given to illustrate practical applications. A classification 
of operator stable Levy processes in two dimensions is provided according to their expo- 
nents and symmetry groups. We conclude with some remarks and extensions to general 
operator self-similar processes. 



1. Introduction 

Self-similar processes form an important and useful class, favored in practical applica- 
tions for their nice scaling properties, see for example the recent books of Embrechts and 
Maejima [H] and Sheluhin et al. [S]. In finance, self-similar processes such as fractional 
Brownian motion and stable Levy motion are used to model prices (or log returns). Recall 
that a stochastic process X = {X(t)} t >o taking values in M. d is self-similar if 

(1.1) {X(ct)} t > ^ {cPX(t)} t > 

at every scale c > 0. Here = indicates equality of finite dimensional distributions, and 
we assume X is stochastically continuous with X(0) = 0. The parameter /3 > is often 
called the Hurst index [14]. Operator self-similar processes allow the scaling factor (Hurst 
index) to vary with the coordinate. Therefore, a process X as above is said to be operator 
self-similar (o.s.s.) if there exists a linear operator B £ GL(M rf ) such that 

(1.2) {X{ct)} t > Q H {c B X(t)} t > Q , 

for all c > 0, where the matrix power c B := exp(i?logc). The linear operator B in (11.21) is 
called an exponent of the operator self-similar process X. If B = (31 for some (3 > 0, then X 
is self-similar. If B is diagonal, then the marginals of X are self-similar, and the Hurst index 
can vary with the coordinate. This is important in modeling many real world phenomena. 
Rachev and Mittnik [35] show that the scaling index will vary between elements of a 
portfolio containing different stocks (see also Meerschaert and Schemer |30j). In ground 
water hydrology, Benson et al. [U [5] use operator self-similar processes to characterize 
spreading plumes of pollution particles. Because the structure of the intervening porous 
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medium is not isotropic, the scaling properties vary with direction, see Meerschaert et al. 
[28] and Zhang et al. [45]. In tick- by-tick analysis of financial data, it is useful to consider 
the waiting time between trades and the resulting price change as a two dimensional 
random vector. Meerschaert and Scalas [32] show that different indices apply to price 
jumps and waiting times, see also Scalas et al. [43]. Results and further references on 
o.s.s. processes can be found in [TTJ Chapter 9] and [27j Chapter 11], see also the pioneering 
work of Hudson and Mason [13] . 

This paper focuses on operator self-similar Levy processes. They belong to the class 
of operator stable processes, which is reviewed in Section [2] Operator stable processes 
admit parametrization by their operator exponents and spectral measures, which is a 
starting point for their analysis. Section [3] presents a method for simulating the sample 
paths of operator stable Levy processes, consisting of a shot noise representation for the 
large jumps, and a Gaussian approximation of the small jumps. Theorem 13.11 justifies 
this method and provides a bound on the error resulting from replacing the small jump 
part by a multidimensional Brownian motion. Section |4] presents a number of examples 
to illustrate the method, including some practical applications. The uniqueness of the 
exponents, which is a critical issue in modeling, is determined by symmetries of such 
processes (see Sections [2] and El). Section [5] provides a classification of operator stable 
Levy processes in two dimensions according to their exponents and symmetry groups. 
There we also give a complete, explicit description of the possible symmetries in terms 
of the exponent and the spectral measure. Using these results, it is possible to construct 
an operator self-similar process with any given exponent and any admissible symmetry 
group. Finally, in Section [6] we provide some concluding remarks and extensions to general 
operator self-similar processes. 



2. Operator stable processes 

We say that a Levy process X = {X(t)}t>o taking values in R d is operator stable with 
exponent B G GL(M <: ') if for every t > there exists a vector b(t) £ M. d such that 

(2.1) X(t) = t B X{l)+b{t) 

where = means equal in distribution. We say that X is strictly operator stable when 
b(t) = for all t > 0. A Levy process is operator self-similar if and only if it is strictly 
operator stable, in which case the exponents coincide [13j Theorem 7]. In general, if X is 
operator stable and 1 is not an eigenvalue of the exponent B, then there exists a vector 
a such that {X(t) — at} t >o is strictly operator stable; a complete description of strictly 
operator stable processes is given by Sato [41]. Henceforth we will always assume that the 
infinitely divisible distribution [i = C(X(1)) is full dimensional, i.e., not supported on a 
lower dimensional hyperplane. The distributional properties of [i determine those of X. 
Indeed, two Levy processes X and Y have the same finite dimensional distributions if and 
only if X(l) and Y(l) are identically distributed. 

A comprehensive introduction to operator stable laws can be found in the monographs 
[16] and [27]. Since an operator stable law is infinitely divisible, its characteristic function 
can be expressed in terms of the Levy representation (see, e.g., [27l Theorem 3.1.11]). The 
necessary and sufficient condition for a d x d matrix B to be an exponent of a full operator 
stable law is that all the roots of the minimal polynomial of B have real parts greater than 
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or equal to 1/2, and all the roots with real part equal to 1/2 are simple, see [16j Theorem 
4.6.12]. Furthermore, we can decompose the space M. d into a direct sum of two -B-invariant 
subspaces M rf = V\ © Vi and write ji = \i\ * fj,2, where \i\ is a normal law supported on V\, 
and H2 is an infinitely divisible law supported on V2 having no normal component. The 
restriction of B to V\ has all eigenvalues with real part equal to 1/2, and the restriction of 
B to V2 has all eigenvalues with real part strictly greater than 1/2. 

The exponent B in (|1.2|) need not be unique. The set of all possible exponents of /i is 
given by [TBJ Theorem 4.6.7]: 

(2.2) £(fi) = B + TS(p) 
where B E £ (/J,) is arbitrary. Here 

(2.3) S(n) := {A e GL(M d ) : A\i = \i * 5 X for some x £ R d } 

is the symmetry group of the probability measure fi, and TS(fi) is the tangent space of 
S(fi) at the identity The tangent space consists of all tangent vectors x'(0) where x(t) is 
a smooth curve on S(fi) with x(0) = I. It is a linear space closed under the Lie bracket 
[A, B] = AB — BA. If S(fi) is finite, then TS(fi) = {0}, and this is the only case when the 
exponent in (|2.ip is unique. If A £ S([i) and B is an exponent of //, then so is A~ X BA. 
When the exponent is unique, we must have AB = BA, so B commutes with S(fi). The 
use of commuting exponents simplifies the analysis of £([J.). Every operator stable law fj, 
has an exponent B c that commutes with 5(/i), see [TBI Theorem 4.7.1]. If /i is operator 
stable with 5(X) = Od, the orthogonal group on M. d , then B c = (51 for some (3 > is the 
only commuting exponent, and \i is multivariable stable with index 1/(3. Since T(Od) = Qd 
is the linear space of skew symmetric matrices, we get from (12. 2h that 

(2.4) £(fi) = (31 + Q d . 

Recall that a matrix Q is skew-symmetric if Q T = —Q, where Q T is the transpose of Q. 

If S(/i) is an arbitrary compact subgroup of GL(R d ), then by a classical result of algebra 
(see, e.g., [6j Theorem 5]) there exists a symmetric positive-definite matrix W and a 
compact subgroup Q of the orthogonal group Od such that 

(2.5) S{n) = w^gw. 
Then ([23]) becomes 

(2.6) £{n) = B + W^HW, 

where TL is the tangent space of Q. 

Theorem 2 in [25] implies that a compact subgroup Q of GL(]R a! ) can be a symmetry 
group of a full dimensional probability distribution on M. d if and only if it is maximal, 
meaning that Q cannot be strictly contained in any other subgroup that has the same 
orbits. For example, the special orthogonal group 0j~ is not maximal because 0\x = OdX 
for every x £ and 0j~ is a proper subgroup of Od- Consequently, cannot be the 
symmetry group of any full dimensional probability measure on M. d . Actually Theorem 2 
in [25] characterizes the strict symmetry group of /i defined by 

(2.7) 5o(^) := {A G GL(R d ) : A\i = n}. 

However, Theorem 5 in Billingsley [6] implies that 5(/i) = Sq(h * S a ) for some a G M. d . 
Hence S(n) must be maximal as well. 
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In this work we will assume that the operator stable law \i has no Gaussian component, 
so that all the roots of the minimal polynomial of B have real parts greater than 1/2. For 
a given exponent B, consider a norm || • \\b on M. d satisfying the following conditions 

(i) for each x G M d , x ^ 0, the map t \— > \\t B x\\B is strictly increasing in t > 0, 

(ii) the map (t, x) i— > t B x from (0, oo) x Sb onto M. d \ {0} is a homeomorphism, 

where Sb = {x G M d : \\x\\b = 1} is the unit sphere with respect to \\-\\b- There are many 
ways of constructing such norms. For example, Jurek and Mason [TU Proposition 4.3.4] 
propose 

(2.8) \\x\\ B = (J \\s B x\\ p s- 1 ds 

where 1 < p < oo and || • || is any norm on M. d . Meerschaert and Schemer [27j Remark 6.1.6] 
observe that if the matrix B is in the Jordan form, then the Euclidean norm satisfies (i)- 
(ii). Moreover, in this case the function 1 1— > \\t B x\\ is regularly varying. Under conditions 
(i)-(ii) we have the following polar decomposition of the Levy measure of an operator stable 
law 

r rco 

(2.9) v(E)= / l E (s B u)s- 2 ds\(du), E€B(R d ), 

Js B Jo 

where A is a finite Borel measure on Sb called the spectral measure of [i. The spectral 
measure is given by 

(2.10) \{F) = v{{x : x = t B u, for some (t, u) G [1, oo) x F}), F G B{S B ) 

and then it follows from (|2.9p and (|2.10p that the spectral measure A is uniquely determined 
for a given Levy measure v, exponent B, and norm The choice of || • \\b is a matter of 

convenience. For example, if B is in Jordan form, then the Euclidean norm || • || is a natural 
choice for || • \\b- Since fi is full, the smallest linear space supporting the Levy measure v 
is M. d (23, Proposition 3.1.20]. Moreover, we have a relation between the symmetries of /i 
and the strict symmetries of v 

(2.11) S{n) = S (v) := {A G GL(R d ) : Au = u}, 

which is valid for any infinitely divisible distribution without Gaussian part. 

3. Accelerated series representation 

Let X = {X(t)} t >o be a proper operator stable Levy process with exponent B, no 
Gaussian component, and characteristic function in the Levy-Khintchine form 

(3.1) logEe^ 1 )) =i(y,x )+ [ (e^ - 1 - i(y, *>1 {W < 1} ) u{dx). 

In this section, we present a practical method for simulating sample paths of this process. 
Our method is based on a series representation [38] in which the small jumps are approxi- 
mated by a Brownian motion [9]. The Gaussian approximation of small jumps accelerates 
the convergence of the series representation, allowing a fast and accurate simulation of 
sample paths. Assume that the Levy measure v is given by (|2.9|) . where Sb is the unit 
sphere with respect to a norm || • \\b satisfying conditions (i)-(ii) of the previous section 
and A is a finite measure on Sb- Our approach to simulation of X is based on a series 
expansion and a Gaussian approximation of the remainder of such series. Such a series 
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expansion falls into a general category of shot noise representations and is a consequence 
of the polar decomposition (|2.9j) . see remark following [381 Corollary 4.4]. Namely, for any 
fixed T > 0, 

(3.2) x(t) = x + ^|l ( o, t ](r J )(^-y) vj-^, te[0,T], 

where {r,} is an iid sequence of uniform on [0,T] random variables, {r^} forms a Poisson 
point process on (0,oo) with the Lebesgue intensity measure, {vj} is an iid sequence on 
Sb with the common distribution X/X(Sb), and 

(3.3) Cj = / / xa r (dx) dr, 

Jj-l J\\x\\<l 

where 

-B \ 



(3.4) Mj1) = p^_J vl€A j 

(see [39l Eq. (5.6)]). The random sequences {Tj}, {r^}, and {vj} are independent. The 
series (13. 2ft converges pathwise uniformly on [0, T] with probability one, see [3H Theorem 
5.1]. Fix e G (0, 1] and define N e = {N e (i)h 6[0 ,T] h Y 

It is elementary to check that N e is a compound Poisson process with characteristic function 
Eexpi{y,N e (t)) =exp^t J £ (e^ y ' sB ^ - l)s~ 2 dsX(du)^ . 

To see this: Observe that the number of terms M e in the sum (|3.5p is Poisson with 
mean 6 t = T\{Sb)/^\ condition on M e = n in the characteristic function, noting that 
(ri/6 e , . . . , Y n /6 e ) is equal in distribution to the vector of order statistics from n IID stan- 
dard uniform random variables; permute the order statistics; and rewrite the characteristic 
function as an integral. Thus N e has the Levy measure 

v e (A)= / l A (s B u)s- 2 ds\(du). 

The remainder 

(3.6) R e (t) = X(t) - N%t), 

is a Levy process independent of N e and -R e (l) has Levy measure v t of bounded support 
given by 

(3.7) v e (A)= [ [ l A (s B u)s- 2 dsX(du). 

Js B Jo 

Therefore, all moments of i? e (l) are finite. A straightforward computation shows that 

(3.8) a t :=Ei? e (l) = x + / xu e (dx)- xu e (dx). 

J\\x\\>l J\\x\\<l 
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Then we have 

X(t) = ta e + N e (t) + {R e (t) - E[R e (t)}}. 
In our main theorem we will show that under certain matrix scaling R e (t) — E[i? e (t)] 
converges to a standard Brownian motion in M. d . Hence any operator stable Levy process 
can be faithfully approximated by the sum of two independent component processes, a 
compound Poisson and a Brownian motion with drift. To this end we will use Theorem 
3.1 in [9]. A simple computation (see [9j Eq. (2.3)]) shows that the covariance matrix S e 
of R e (l) is given by 



S e = E 

(3.9) 



(R e (l)-E[R t (l)])(R t (l)-K[R e (l)})~ 



where A is given by 



(A)(A)'r^A((iu)= / s^A(s^) 1 s-'ds, 
s B Jo Jo 



(3.10) A= / uu 1 X(du). 

Js B 

We observe the following scaling 

(3.11) S e = e- 1 [\er) B A((er) B ) T r~ 2 dr = e~ 1 e B J: 1 (e B ) T . 

Jo 

Theorem 3.1 in [9] assumes that £ e is nonsingular for all e > 0. Since we consider the 
spectral measure together with the exponent B as primary parameters of an operator stable 
law, it is natural to state the nonsingularity condition in terms of these characteristics. 
Let line(suppA) denote the smallest .B-invariant subspace of M. d containing the support of 
A. If v is as in (|2.9p . then the support of u is not contained in a proper subspace of W 1 if 
and only if 

(3.12) lin B (supp A) =R d 

cf. [16], Corollary 4.3.5. In particular, (|3.12|) holds when A is not concentrated on a proper 
subspace of M. d . 

As we have stated in Section [U since X does not have Gaussian component, 

(3.13) K ■= min{6i, ...,b d }>-, 

where b±, . . . , bd are the real parts of the eigenvalues of B. This quantity does not depend 
on a choice of B because the real parts of eigenvalues of all exponents of X are the same 
(see [271 Corollary 7.2.12]). 

Theorem 3.1. Let X be an operator stable Levy process with exponent B and let the Levy 
measure of X(l) be given by (|2.9p such that (|3.12l) holds. Fix T > and let N e be as in 
(|3.5p . W be a standard Brownian motion in M. d independent o/N e , and a e = {a e i}t>o be 
a drift determined by (|3.8p . Define 



(3.14) A e = e- l ' 2 e B T^ 2 



where Si is given by (|3.9I) with e = 1. 

Then, for every e £ (0, 1] there exists a cddldg process Y e such that on [0, T] 

(3.15) X = a e + A e W + N e + Y e 
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in the sense of equality of finite dimensional distributions and such that for every 5 > 

(3.16) e l/2-6.+* gup ||y e ( t )|| flse ^ 

te[o,T] 

where is given by (|3.13p . 

Proof. First we will prove that Ei is nonsingular. Let v\ be the Levy measure (13, 7p with 
e = 1 and let 

L = lin(supp v\) 

be the closed linear space spanned by suppz^i. By [H Lemma 2.1] it suffices to show that 
lin(supp^i) = R . Following [HI Corollary 4.3.5] we have 

suppz/i = {x : x = s B u, < s < 1, u G suppA}. 

We will show that L is U-invariant. To this end it is enough to show that if x = s B u G 
suppz^i, for some < s < 1 and u G suppA, then Bs B u G L. For any £ (0,1), 
(9s) B u G supp i/i so that 

„ o , ((9s) B u - s B u T 

Bs B u = lim ^— i- G L. 

0/1 log 

Since L is closed and -B-invariant and contains the support of A, L = M. d by (|3.12|) . Thus 
Si is nonsingular. 

Theorem 2.2 in [9] shows that the asymptotic normality of R e (t) —K[R e (t)] holds if and 
only if for every At > we have 

(3.17) lim/ (E- 1 x,x)u e (dx) = 0. 
Using (13. lip we have 

(Z: l s B u, s B u) = 6(( e - B ) T Sr 1 e- B S B n, s B u) 
= e(^ 1 e- B s B u,e- B s B u) 
= e(^ 1 (s/e) B u, (s/e) B u) 

Note that in general (Ax,x) < \\A\\\\x\\ 2 < C\\ A\\ \\x\\g (for some constant C > 0, since all 
norms on M. d are equivalent). Then, since t i— > is strictly increasing and t B x = x 

when t = 1, the above bound shows that 

(3.18) {E- 1 s B u,s B u) < CellS^IIIKs/e) 5 ?/!!! < CeH^^ 1 1|, 

whenever < s < e < 1 and u G Sb- Since Si is invertible we know that c\ = C\\Y,^ 1 \\ G 
(0, oo). Then, for every k > and e G (0, 1) we have 



(E e l x, x) v e {dx) 

<E7 1 a B u, s B n) s~ 2 ds\(du) 

{(s,ti)e(0,e]xS B : (S7 1 s fl u,s B ii)>K} 

= 

when e < c{ l K. Indeed, in view of (I3.18P the region of integration is empty for c\e < k. 
Therefore, (|3.17l) trivially holds. 
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Applying [9j Theorem 3.1] we get (13.15P and that 

(3.19) sup H^r 1 !^*)!! ase^O. 

te[o,T] 

It remains to show (13. 16p . If ||Ei|| = C2 then HE/ || = Since every eigenvalue of —2? 
has real part less than or equal to — [271 Proposition 2.2.11 (d)] implies that for any 
5 > 0, for some C3 > 0, we have ||t~ B x|| < c 3 ^ fe * +<5 ||x|| for all t > 1 and all x G M d . Then 
< ess 6 *" 15 for all a < 1. Then for all < e < 1 we have 

\\M ^e-v^n^iiiisj/ 2 !! < ce -V2-5+6,. 

where c = 03^/02. Therefore, 

||y e (t)|| < UMA^Um^ce^-^llA^YMl 
which together with (13, 19ft yields (|3.16l) . The proof is complete. □ 



4. Simulation 

The main goal of this paper is to provide a practical method for simulating the sample 
paths of an operator stable Levy process X. Theorem 13.11 decomposes X into the drift 
a e t, the large jumps N e (t), and a Gaussian approximation of the small jumps, with a 
remainder term whose supremum converges to zero in probability at a polynomial rate as 
the number of large jumps increases (or, equivalently, as the size of the remaining jumps 
tends to zero). In this section, we will demonstrate the practical application of Theorem 
13.11 and illustrate the resulting sample paths. 

Theorem [3J] justifies the use of the process 

(4.1) Z e (t) := a e t + A e W(t) + N"(t), 

with A e given by (13. 14ft and W(t) a standard Brownian motion, to simulate sample paths of 
the operator stable process {X(t)} te [ ^T] specified by (|3.ip and (|2.9p . Formula (I3.16P shows 
that the approximation converges faster when the real parts of the eigenvalues of B are 
uniformly larger. Remark 7.2.10 in [27] shows that the real parts of the eigenvalues of the 
exponent B govern the tails of the operator stable process X(t), and 6* = min{6i, . . . , bd} > 
1/2 determines the lightest tail, in the sense that E,\(X(t),u)\ p diverges for all p > 1/6* 
and all u 7^ 0. Hence the convergence is faster when X has a heavier tail. 

The process Z e (t) in (14. ip approximates the operator stable process {X(t)} te \o t T] by 
discarding small jumps, replacing their sum by an appropriate Brownian motion with 
drift. The discarded random jumps are all of the form r B v where v G Sb and r < e. If B 
has no nilpotent part then ||r B ?;||£ < e b * . Hence in order to retain all jumps larger than 
m it suffices to take e = m 1 ^* , and then the number of jumps simulated will be Poisson 
with mean m~ l l b *T\(SB)- If there is a nilpotent part, the bound involves additional loge 
terms. 

In general, an operator stable process can be decomposed into two independent com- 
ponent processes, one Gaussian and another having no Gaussian component. The two 
components are supported on subspaces of M d whose intersection is trivial. In practical 
applications, Theorem 13. II is applied to the nonnormal component. In the case where X(t) 
has both a normal and a nonnormal component, the resulting approximation combines a 
full dimensional Brownian motion with drift, and a Poissonian component restricted to 
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the nonnormal subspace. For the remainder of this section, we will focus on simulating 
operator stable laws on R 2 having no normal component. 

In practical applications, it is advantageous to produce a simulated process whose mean 
(if the mean exists) equals that of the operator stable process X(t). If every eigenvalue of 
the exponent B has real part b < 1, then the mean exists, by [27j Theorem 8.2.14]. If any 
eigenvalue has real part b > 1 then the mean is undefined. In the former case, one can 
choose a e so that the right-hand side in (|4.ip has mean zero. Recall that the number of 
terms M e in the sum (|3.5p defining N e (t) is Poisson with mean 9 e = T\(Sb)/c, and that 
conditional on M e = n, (Ti/9 e , . . . ,T n /9 e ) is equal in distribution to the vector of order 
statistics from n IID standard uniform random variables. Condition to get K[N e (t)\M e = 
n] = n(t/T)K[(eU)~ B ]E[v] where U is standard uniform and v has distribution A/A(«S#). 
Removing the condition and simplifying shows that 

(4.2) E[N e (t)] = t\{S B )e B - I E[U- B }E[v}. 

Since E[W(t)] = we can set a e t = — E[iV e (t)] to get mean zero. Note that for such B we 
have ||e s_/ o;|| — ► oo for all x / by [27l Theorem 2.2.4], so that ||a e || — ► oo as e — ► 0. This 
reflects the fact that, in the finite mean case, the infinite series (|3.2p does not converge 
without centering. Finally we note that, if K[v] = 0, then no centering is necessary. 

In this section, we assume a fixed coordinate system on K 2 with the standard coordinate 
vectors e\ = [1,0] T and e 2 = [0, 1] T , and we write X(t) = X\(t)e\ + X2(t)e2- Recall that 
a strictly operator stable process satisfies the scaling relationship 

(4.3) X(t) = t B X(l) 

for all t > 0. All plots in this section use T = 1 and e = 0.001, and we show the simulated 
processes at the time points t = nAt for < t < T with At = 0.001. Unless otherwise 
noted, we use the standard Euclidean norm. 

Example 4.1. Equation (14. ip was used to simulate an operator stable process X(t) whose 
exponent is diagonal 

" 1/1.8 
1/1.5 

so that Bei = biei with b\ = 1/1.8 and 62 = 1/1.5. Since the exponent is already in Jordan 
form, we can take to be the usual Euclidean norm, so that Sb is the unit circle. 

We choose the spectral measure A to place equal masses of 1/4 at the four points ±ei 
and ±e2- Then K[v] = in (14. 2h so that no centering is needed, as the simulated process 
has mean zero without any centering. Then A = diag(l/2, 1/2), Si = diag(9/2, 3/2), and 
A e = diag(3\/5v/l0/10, v^/10). It is easy to see from the definition t B = I + Blogt + 
(Blogt) 2 /2l H that t B = diag(t fel , t b2 ). From the scaling relation ()43|) it follows that 

X(t) = i 6 *AQ(l). 

Hence the coordinate marginals are (strictly) stable with index ot\ = l/b± = 1.8 and 
a.2 = 1.5, respectively. The top right panel in Figure 1 shows a typical sample path of the 
process, an irregular meandering curve punctuated by occasional large jumps. The top left 
panel shows the corresponding shot noise part N e (t) before the Gaussian approximation 
of the small jumps is added. Since the spectral measure is concentrated on the coordinate 
axes, the large jumps apparent in the sample path of Figure 1 are all either horizontal or 
vertical. Pruitt and Taylor [34] showed that the Hausdorff dimension of the sample path is 
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Figure 1 . Simulated operator stable process for Example 14.11 with inde- 
pendent symmetric stable marginals. Top left panel shows the sample path 
of the shot noise process M e (i), and top right panel shows the correspond- 
ing operator stable process X(t). Bottom left panel shows the marginal 
process Xi(t), and bottom right panel shows X^if)- 



maxjai, a^} = 1-8 with probability one. Since the spectral measure is concentrated on the 
coordinate axes, Lemma 2.3 in Meerschaert and Schemer [26] shows that the coordinates 
X\{t) and X2(t) are independent stable processes. The bottom panels in Figure 1 graph 
each marginal process. Note that the large jumps occur at different times, reflecting the 
independence of the marginals. Blumenthal and Getoor [7j showed that the graph of the 
stable process Xi(t) has Hausdorff dimension 2 — aj. 

Example 4.2. The same exponent B is used as in Example 14.11 but now we take the 
spectral measure A(ej) = 1/2. The matrices A and A e turn out to be the same as Example 
14.11 The marginals Xi(t) are still stable with index a\ = 1.8 and 02 = 1.5, but they 
are no longer symmetric, and we center to zero expectation. From (|4.2|) we get a e = 
[45^10/4, 15] T to compensate the shot noise portion to mean zero. Figure 2 shows a 
typical sample path and component graphs for this process. Since the spectral measure is 
concentrated on the positive coordinate axes, the large jumps apparent in the component 
graphs are all positive. 
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Figure 2. Simulated operator stable process for Example 14.21 with inde- 
pendent skewed stable marginals. Top left panel shows the sample path of 
the shot noise process M e (t), and top right panel shows the corresponding 
operator stable process X(t). Bottom panels show the marginal processes 
Xi(t) and X 2 (t). 

Example 4.3. The same exponent B is used as in Example 14. 11 but now we take the spec- 
tral measure A to be uniformly distributed on the unit circle: set v = (x 2 + y 2 )" 1 / 2 ^, y] T 
where x, y are independent standard normal. The matrices A and A e turn out to be the 
same as Example 14.11 Since K[v] = 0, no centering is needed. The marginals Xi(t) are 
symmetric stable with index ot\ = 1.8 and «2 = 1.5, but they are no longer independent. 
The top panels in Figure 3 show a typical sample path of the process. Since the spectral 
measure is uniform, the large jumps apparent in the sample path take a random orienta- 
tion. Theorem 3.2 in Meerschaert and Xiao [H] shows that the sample path is a random 
fractal, a set whose Hausdorff and packing dimension are both equal to 1.8 with probability 
one. The bottom panels in Figure 3 show the graphs of each marginal process. Note that 
the large jumps in both marginals are simultaneous, reflecting the dependence. 

Example 4.4. Figure 2 of Zhang et al. [15] represents a model of contaminant transport in 
fractured rock. Pollution particles travel along fractures in the rock, which form at specific 
angles due to the geological structure of the rock matrix. An operator stable process X{t) 
represents the path of a pollution particle, with independent skewed stable components in 
the fracture directions. The skewness derives from the fact that particles jump forward 
(downstream) when mobilized by water that flows through the fractured rock. The two 
components of X(t) are skewed stable with index a = 1.3 on the line with angle B\ = 30° 
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Figure 3. Simulated operator stable process for Example 14. 3 1 with depen- 
dent stable marginals. Top left panel shows the sample path of the shot 
noise process, and top right panel shows the corresponding operator stable 
process. Bottom panels show the marginal processes. 




Figure 4. Simulated operator stable sample path for Example I4.41 with 
independent skewed stable components along nonstandard coordinate axes. 
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measured from the positive e± axes as usual, and index 1.7 on the line with angle 62 = —35°. 
The two stable laws are independent. The e\ axis represents the overall direction of flow, 
caused by a differential in hydraulic head (pressure caused by water depth). The exponent 
B has one eigenvalue b\ = 1/1.3 with associated eigenvector V\ = Ro 1 e\ = [.865, .500] T , 
and another eigenvalue 62 = 1/1-7 with associated eigenvector V2 = Re 2 e i = [-820, — .572] T . 
The spectral measure is specified as X(v\) = 0.4 and X(v2) = 0.6, representing the relative 
fraction of jumps along each fracture direction. In order to compute the matrix power t B 
a change of basis is useful. Define the matrix P according to Pej = V{ so that 



P 



.865 .820 
.500 -.572 



and D = P 1 BP = diag(&i, 62) is a diagonal matrix. Then the exponent 



B = PDP~ 



From (|3~T0]) we get 



A 



.703 
-.109 



.688 .142 
.057 .671 



-.109 
.297 



Since t D = diag(i fel , t b2 ) we can compute t B = Pt D P 1 and integrate in (13. 9ft to get the 
Gaussian covariance matrix S e whose symmetric square root is given by 

0.723 -0.416 
-0.416 0.407 



To compute the square root, we decompose S e = QEQ 1 where E = diag(ci,C2), Cj are 
the eigenvalues of E e , and the columns of Q are the corresponding eigenvectors, so that 



A e = QE 1 / 2 Q~ 1 where E 1 ' 2 = diag(c 



1/2 



From (|4.2|) we get a e 



to compensate the shot noise portion to mean zero. Note that B 



II; 



= [27.9,-10.1] T 
biUi where u\ = 



[.572, .820] 1 and U2 = [.500, — .865] T are the dual basis vectors. Then each projection 
(X(t),Ui) is (strictly) stable with index oti = 1/fej, since 



(X(t), Ui ) i (t B X(l), Ul ) = (X(l),t BT Ul ) = (X(l),t b > Ui )=t b >(X(l), Ui ). 

Hence .572X x (t) + .820X 2 (t) is stable with index a x = 1.3 and .500Xi(t) - .865X 2 (t) 
is stable with index «2 = 1-7. Lemma 2.3 in [26] shows that these two skewed stable 
marginals of X(t) are independent, since the spectral measure is concentrated on the 
eigenvector coordinate axes (x,Vi) = 0. Figure 4 shows a typical sample path, along with 
the coordinate marginals. Note that the large jumps lie in the i>, directions. The mean zero 
operator stable process X(t) represents particle location in a moving coordinate system, 
with origin at the center of mass. Hence Figure 4 illustrates the dispersion of a typical 
pollution particle away from the center of mass of the contaminant plume. Dispersion is 
the spreading of particles due to variations in velocity, and it is the main cause of plume 
spreading in ground water hydrology. 



Example 4.5. We simulate an operator stable process X(t) whose exponent has a nilpo- 
tent part 

" 1/1.5 
1 1/1.5 



B 



1 1 
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We choose the spectral measure A to place equal masses of 1/4 at the four points ±ei and 
±e2. Then E[i>] = in (14. 2p so that no centering is needed. Here A = diag(l/2, 1/2), 



3/2 
-9/2 



and 



0.146 
-0.359 

Note that t B = t b t N where b = 1/1.5 and 



-9/2 
57/2 

-0.359 
4.009 



t N = 



1 

lost 1 



From (14.31) it follows that the second marginal X%{t) is symmetric stable with index a = 
1/6 = 1.5. The first marginal is not stable, but it lies in the domain of attraction of a 
symmetric stable with index a = 1.5, see pH Theorem 2]. Figure 5 shows a typical sample 
path of the process. The large jumps apparent in the sample path of Figure 5 are all of 
the form t B v where v = ±ej and t > 0. Hence they are either vertical, or they lie on the 
curved orbits ±t B e\. Theorem 3.2 in [31] shows that the sample path is almost surely a 
random fractal with dimension 1.5. Lemma 2.3 in [26] shows that the coordinates X±(t) 
and X2(t) are not independent. 

Example 4.6. We simulate an operator stable process X(t) whose exponent 

~ 1/1.5 1 
-1 1/1.5 
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Figure 6. Simulated operator stable sample path for Example I4.61 whose 
exponent has complex eigenvalues. 



has complex eigenvalues b ± i with b = 1/1.5. We choose the spectral measure A to place 
equal masses of 1/4 at the four points ±ei and ±e2, so that K[v] =0 in (14. 2ft and no 
centering is needed. Here A = diag(l/2, 1/2), and 

0.387 0.0 
0.0 0.387 

In this case, t B = t b R 9 ^ with 9(t) = lni, since we can write B = bl + K where the matrix 
exponential exp(si^) = R s . The coordinate marginals X\{t) and X2(t) are not stable, but 
they are both semistable with index a = 1/6 = 1.5, see [2H Theorem 2]. Lemma 2.3 in 
[26] shows they are not independent. Figure 6 shows a typical sample path of the process. 
The large random jumps are of the form t B v where v = ±ej, so that the angle varies along 
with the length of the jump. The sample path is a fractal with dimension 1.5, see [3"T1 
Theorem 3.2]. 

Example 4.7. Figure 1 in Zhang et al. [45] presents an operator stable model X(t) with 
diagonal exponent 

" 1/1.5 
1/1.9 

and spectral measure A that places masses of 0.3 at e±, 0.2 at ±6°, 0.1 at ±12°, and 0.05 
at ±18° on the unit sphere in the standard Euclidean norm. Large jumps are along the 
positive x-axis, or along the orbits t B u where u is a unit vector at ±6°, ±12°, or ±18°, 
representing displacements of a pollutant particle in an underground aquifer with a mean 
flow in the positive x direction, but some dispersion due to the intervening porous medium. 
The average plume velocity is v = [10, 0] T so that E[X(t)] = tv. Figure 7 depicts the path 



B 
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Figure 7. Simulated operator stable process for Example 14.71 modeling a 
pollution particle moving through underground water in a heterogeneous 
porous medium consisting of sand, gravel, and clay. Left panel depicts the 
sample path of a moving particle. Right panel shows the coordinate 
marginals. 



of a typical particle. Here 



and 



A 



0.977 
0.0 



0.0 
0.0226 



0.541 0.0 
0.0 0.546 

From (|4.2I) we compute a e = [29.7, 0] T and, in the simulation code, we first center to mean 
zero, and then add the mean velocity. 



Example 4.8. This example follows the transport model number 22 for contaminant 
transport in complex fracture networks from Reeves et al. [37]. The exponent B has 
eigenvectors v x = [y/2/2, V2/2] T and v 2 = [y/2/2, -y/2/2] T at +45° and -45° on the unit 



circle with eigenvectors b± = 1/1.1 and b 2 
D = P- 1 BP = diag(6i, b 2 ) so that 

B = PDP^ 1 = 



1/1.2 respectively. Writing Pet = V{ we get 



115/132 5/132 
5/132 115/132 



Since t D = diag(i bl , t b2 ) we also have 

t B = pt Dp-l = 



(t bl +t b2 )/2 (t bl -t b2 )/2 
{t bl -t b2 )/2 (t bl +t b2 )/2 



The spectral measure has weights 0.4 and ±45° and 0.2 at e 2 - The Levy measure is 
concentrated on the two straight line orbits {t B Vi : t > 0} and on the curved orbit {t B e 2 : 
t > 0}. Marginals (X(t),Vi) are stable with index ot\ = 1.1 and a 2 = 1.2 respectively, but 
they are not independent, since the spectral measure is not concentrated on the eigenvector 
axes. The first marginal process (X(t),vi) is positively skewed, since the projection of the 
Levy measure onto the first eigenvector coordinate places all mass on the positive half line. 
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Figure 8. Simulated operator stable sample path for Example 14.81 mod- 
eling the motion of a pollution particle moving through water in fractured 
rock. The mean zero sample path represents deviation from the plume 
center of mass. 



The second marginal (X(t),V2) is the sum of two independent stable processes, one with 
positive skewness resulting from the V2 orbit, and one with negative skewness resulting from 
the projection of the e 2 orbit onto the negative t>2 axis. As in Example 14.41 we compute 
A = diag(0.4, 0.6) and 

0.0603 -0.0204 



-0.0203 0.0723 



From (14. 2} we compute a e = [11.35, 4. 43] T to correct the shot noise process to mean zero. 
Figure 8 shows a typical sample path. In this case, the sample path represents the growing 
deviation of a typical pollution particle from the plume center of mass. 



Example 4.9. This example is an operator stable process whose exponent 

_ r 1/1.8 1/2 

[ 1/1.5 ' 

The matrix B has eigenvalue-eigenvector pairs Bvi = bivi with b\ = 1/1.8, v\ = e\, 
62 = 1/1-5, and vi = (9/2)ei + ei. As in Example 14.41 we compute 

, B _ [ *Vi-8 -(9/2)^/1-8 + (9/2)^/1-5 - 

[ fVi-5 

We choose the norm ((23]) with p = 2. Compute \\x\\ 2 B = (9/10)xf - (81/110)xix 2 + 
(903/880)^2 so that the unit sphere Sb is an ellipse, whose major axis is rotated ap- 
proximately 50° counterclockwise from the e\ direction. Figure 9 shows level sets of this 
norm. The spectral measure A places equal masses of 1/4 at each point where the unit 
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Figure 9. Level sets of the norm \\x\\b used in example 14.91 



sphere Sb intersects the coordinate axes: ic^ where c\ 
A = diag(5/9, 440/903), and 



10/9 and c 2 , = 880/903. Here 



5.209 
-0.266 



-0.266 
0.274 



The second coordinate ^(i) is symmetric stable with index a?, = 1.5, and the projection 
onto the remaining eigenvector X±(t) — (9/2)^2 (t) is stable with index ot\ = 1.8. These two 
stable marginals of X(t) are not independent, since the spectral measure is not concentrated 
on the eigenvector axes. Figure 10 shows a typical sample path for this process. The large 



jumps of the process are all of the form t v where v = ±ciei or v = ±C2e2, since we 
have concentrated the spectral measure at these points. If v = ±ciei then, since e± is an 
eigenvector of B (and hence of t B ), these jumps will be in the horizontal. The remaining 
jumps lie along the orbits {±t s C2e2 : t > 0}. 



5. Exponents and Symmetries in two dimensions 



Section [4] illustrates the wide range of possibilities represented by operator stable Levy 
processes in M 2 . In this section we will provide a classification of such processes, according 
to the type of exponent, and the symmetry group. Let X = {X(t)}t>o be an operator 
stable Levy process with exponent B, and fx = C(X(1)), as in Section [2j Since the real 
parts of the eigenvalues of B are greater than 1 /2, after a change of coordinates, if needed, 
the exponent assumes one the following Jordan forms 



(5.1) 



B = bl, B 1 



bi 
b 2 



B, 



b 
1 b 
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_ 4 i , 1 _ 2 i , 

0.5 1 0.5 1 

Figure 10. Operator stable sample path for example 14.91 whose exponent 
B is not given in Jordan form. 



where 6,61,62 > 1/2, 61 7^ 62, and c / 0. If B = Bq, then X is a multivariable stable 
process with index a = 1/6, and all maximal compact subgroups of GL(M 2 ) are admissible 
as S(n). A genuine operator stable Levy process is obtained when B = Bi, i = 1, 2, 3. Our 
first question is, what are possible symmetry groups? 

To deal with this question, we need to review some basic facts about subgroups of the 
orthogonal group 02 on M 2 , which can be found, e.g, in [2]. Recall that O2 consists of 
rotations and reflections, 

2 = {R e ,Fg:6e [0,2tt)}, 



where 



Re 



cos 9 
sin 9 



— Sill! 

cos 9 



and Fg 



cos 9 sin 9 
sin 9 — cos t 



Re is a rotation counter-clockwise by 9 and Fg is a reflection through the line of angle 9/2 
passing through the origin. The following rules of composition hold: Rg 1 Rg 2 = Rg 1+ g 2 , 
Fe 1 Fg 2 = Re 1 -e 2 i ReiFe 2 = Fg 1+ g 2 , Fq 2 R 9i = Fq 2 _ 9i . 

The group of rotations 0^ = {^8 '■ @ e [0,27r)} is the only infinite proper compact 
subgroup of 02- There are also only two kinds of finite subgroups of O2 (modulo the 
orthogonal conjugacy, see [21 Ch. VII. 3]): 
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(1) Cyclic group C n = {R k27T /n ■ k = 0, . . . ,n - 1}, n > 1, 

(2) Dihedral group V n = {R k2n / n , F k2n/n : k = 0, . . . ,n - 1}, n>\. 

Notice that d = {I}, C 2 = {I, -I}, V x = {I, F }, and V 2 = {I, F Q , -I, -F }, where 



F 



1 








is the reflection with respect to the x-axis. We will also need T>\ = {I, — Fq}, the group of 
reflection with respect to the y-axis, which is orthogonally conjugate to T>\. 

The next result characterizes the possible symmetries of the distribution of X(t) in the 
truly operator stable case where B = Bi in (|5.ip for some i = 1, 2, 3. In view of (12.11) . the 
symmetry group do not depend on t. Remarkably, once the exponent takes the Jordan 
form, all symmetries must be orthogonal, not just conjugate to an orthogonal matrix. 

Theorem 5.1. Let X = {X(t)}t>o be a full operator stable Levy processes on M 2 with an 
exponent B in the Jordan form (15. ip . and let fi = C(X(1)). Then the following hold. 

(i) If B = B\, then S(fi) is either C\, C 2 , T>\, T>\, or T> 2 . 

(ii) If B = B 2 , then S(fi) is either C n , n > 1, or 2 . 
(hi) If B = i?3, then S(fi) is either C\ or C 2 . 

Proof. Suppose that [i has an exponent B = Bi, i = 1,2,3, and let B c be a commuting 
exponent, see Section [2l If S{pi) is finite, then B. L = B Cl otherwise B c can be different 
from Bi. The symmetries 5(/i) defined in (|2.3|) form a compact subgroup of the centralizer 

C(B C ), 

(5.2) S(fj,) C C(B C ) := {A G GL(M 2 ) : AB C = B C A}. 

First consider finite symmetry groups S(/j,), so that B c = Bi. If i = 1, 

a 
p 

and since S(fi) is finite (and thus compact), 

a 



C{B 1 







a 



1 



Thus 5(/i) is either C\, C 2 , T>\, T>\, or V 2 , as claimed. If i = 2, then 



C(B 2 ) 

Since S(fi) is finite (and thus compact) 

5(m)c 



-0 



a 



a 



so that S(n) = C n , for some n > 1. If i = 3, 

a 

(3 a 



C(B 3 ) -- 

and since 5(/i) is finite, 
Thus S(fi) is either C\ or C 2 , as claimed 



a 2 + (3 2 > 



a 2 + (3 2 = 1 



{[Sa]=M = l}- 
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Now we consider infinite symmetry groups <S(//), so that <S(/x) = W~ 1 2 W for some 
symmetric positive definite matrix W, see (|2.5|) . From (|5.2p . WB c W~~ l commutes with 
every orthogonal transformation. Thus WB c W~ l is a multiple of the identity matrix, 
which yields 

(5.3) B c = (31 

Since T0 2 = Q 2 , 

Bi = B c + W~ l KW = W~ l {fiL + K)W 
for some skew symmetric matrix K, and so 

for some 7 / and <p G [0, 27r). This equation eliminates the cases i = 1 and i = 3 by 
comparing the eigenvalues on the left and right hand side. Thus i = 2 and £?2 = a-R^ f° r 
some ip G (0, tt) U (7r,27r), from which we have 

Comparing the determinants of both sides gives a = 7. Hence 

ify = W^R+W. 

Since the sets of eigenvalues of both sides of this equation must be the same, <p = ip or 
4> = 27r — tp. If (j) = tp then VKi?^, = R^W for -0 G (0,vr) U (vr,27r). A direct verification 
of this equation reveals that W = nR T is a multiple of a rotation. (In fact, W is a scalar 
multiple of the identity, since it is also symmetric and positive definite.) Therefore, 

S(ji) = (kR t )- 1 2 kR t = 2 , 

as claimed. If <\> = 2tt — ip, then 

ify = W~ 1 R 2n ^W = W^FqF^W = W^FqR^FqW 

or 

(F W)R >P = R^(FoW). 
By the same reason as above, one can verify that FqW = kR t is a multiple of rotation. 
Hence W = kF„ t and 

S(fl) = (kF„ t )" 1 02^t = 2 . 
This proves that B = B 2 and S(ji) = 2 provided S{p) = W~ l 2 W. □ 

Remark 5.2. Any operator stable law can be transformed to one in which the exponent 
takes the Jordan form (I5.1|) by a simple change of basis. Theorem I5.ll shows that the 
Jordan basis renders all symmetries orthogonal, and then the Jordan form determines 
which symmetries are possible. 

Operator stable laws are parameterized by their exponents and spectral measures. 
Therefore, it is useful to have their symmetries described in terms of these parameters. 
Recall that O2 denotes the subgroup of rotations of the orthogonal group 2 . 

Theorem 5.3. Let X = {X(t)}t>o be a full operator stable Levy process in M 2 with 
exponent B and no Gaussian component, and let fi = C(X(l)). Suppose that B is 
given in the Jordan form (|5.ip and that the spectral measure A is determined by the 
polar decomposition (|2.10p relative to Sb = S , the Euclidean unit sphere of M 2 . Let 
<Sq(A) = {A G GL(IR 2 ) : AX = A} denote the strict symmetry group of the spectral measure. 
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(a) If B = B 1 , then S(ji) = S (A) n V 2 . 

(b) If B = B2, then either S(fi) = So (A) n = C n for some n > 1, or S(fi) = 
S (X) = O 2 . 

(c) If B = B 3 , then S(p) = S (X) n C 2 . 

Proof. Let v be the Levy measure of [i. Since \x does not have a Gaussian part, we have 

(5.4) = S (u) = {Ae GL(M 2 ) : Av = u} 

as in (|2.11|) . First we will show that if B = Bi, i = 1, 2, 3 and S(/i) is finite, then 

(5.5) S(p) = S (A) n{A£0 2 :AB = BA}. 
Indeed, recall (23) and d^TTO) for S B = S 1 : 

r roc 

v(E) = / / l £ (s jB u)s- 2 dsA(dn), £ G 5(M 2 ), where 
Js 1 Jo 

X(F) = v{{x : x = t B u, for some (t, it) G [1, 00) x F}), F G ^(S 1 ). 

Let ^4 G S(fi), S(fi) being finite. Then A G 2 by Theorem 15.11 and ^4 commutes with B. 
For every F G B{S l ), A~ l F G ^(S 1 ) and 

A(A" 1 F) = z^({x : x = t B A~ l v, for some G [l,oo) x F}) 

= v(A~ l {x : x = t B v, for some (t, v) G [1, 00) x F}) = X(F) 

because <S(/x) = So(u) from |5.4jl . Hence ^4 G So (A). The proof of the opposite inclusion 
in (15.51) uses similar arguments and is omitted. 

Proof of (a). A direct verification shows that B\ commutes with V 2 . Thus by (|5.5p 

<S (A) D%) DSo(A)nP 2 . 

Since S(/.i) C P2 by Theorem 15.11 we get (a). 

Proof of (b). By Theorem 15.11 S(p) = C n for some n > 1, or S(/i) = 02- Suppose that 
<S(//) = C n . Since O2" commutes with B 2 , by (|5.5p we have 

S (A) D SO) D S (A) n 0+ 

Thus S(/i) = S (A) n = C n . 

Suppose S(/i) = 2 . Then Rg G So(^) for every 6 by (15.41) . Since i?g commutes with 
B 2 , Rq G So (A) by the same line of arguments as in the proof of (|5.5p . Hence So (A) D > 
which implies that A is a finite full measure in IR 2 . Then A is a constant multiple of a 
probability measure, so So(A) must be maximal by [25, Theorem 2], and hence So (A) = 2 . 

Proof of (c). It follows from (15. 5|) because C 2 obviously commutes with B3. □ 

Remark 5.4. With the help of Theorem 15.31 it is possible to explicitly construct an op- 
erator stable process with any given exponent Bi for i = 1,2,3 in the Jordan form (|o.ip 
and any admissible symmetry group. For example, let A be concentrated at four points 
(±2 -1 / 2 , ±2 -1 / 2 ). Choosing masses at these points appropriately, any subgroup of T> 2 is 
realized as So(A). By Theorem 15.31 all cases of S(/z) are realized by this example when 
B = B\ and B = B3. When B = B 2 , we only get C\ and C 2 . To get S(/z) = C n , n > 3, we 
take A concentrated at vertices of a regular n-gon inscribed into the unit circle with one 
vertex at (1,0) and equal masses at all the vertices. Then So (A) = V n , so by Theorem 15.31 
S(/i) = V n n O2 = C n - S(/i) = 2 when B = B 2 and A is a uniform measure on S 1 . 
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Figure 1 1 . Support of the Levy measure (thick lines) for Remark 15.51 
showing the effect of the exponent B\ on the symmetry group. The spectral 
measure gives equal weight to three equally spaced points on the unit circle, 
so that «So (A) =T>3. In this case B = B\, we have S{pi) = T>±. 



Remark 5.5. It is interesting to see how much an exponent affects the symmetry. Consider 
a measure A with <So(A) — Dn described in Remark 15.41 with n = 3. Then, by Theorem 
15.31 = T>\ when B = B\, <S(/i) = C3 when B = B2, and <S(/i) = C\ when B = B3. 

Figure 11 illustrates the diagonal case B = B\, in which the Levy measure v in (|2.9p is 
symmetric with respect to reflection about the vertical axis. Here we take 61 = 1/1.8 and 
62 = 1/1-5, but any case with b\ 7^ 62 appears similar. Figure 12 illustrates the complex 
case B = B2, where the Levy measure is symmetric with respect to rotations that are 
a multiple of 27r/3. Figure 13 illustrates the nilpotent case B = B3, and here the Levy 
measure has no nontrivial symmetries. All three cases have the same spectral measure, but 
a different exponent. Hence the spectral measure and the exponent are both important in 
determining the symmetries. 

Remark 5.6. In order to tie the theoretical results of this section back to the concrete 
examples in Section [4j we compute the symmetry group <S(//) for a few interesting cases. 
For Example 14.11 we have S(fx) = <So(A) = T>2 by Theorem 15.31 (a), since the exponent 
B = B\ in (15. ip . and spectral measure A gives equal mass to the four points ±ei,±e2. 
The spectral measure in Example 14.31 is uniform on the unit sphere, so that <So(A) = O2, 
but the symmetry is of the form B = B\ in (|5.ip . so the symmetry group <S(/i) = T>2 by 
Theorem 15.31 (a). The construction in Example 14.51 yields <So(A) = T>2- Then <S(/i) = C2 
since the exponent B = B3 is nilpotent, by Theorem 15.31 (b). In Example 14.61 we also have 
So(A) = ^2, and then <S(/x) = C2 by Theorem 15.31 (c). Example 14.71 has <S(/i) = T>\ since 
the spectral measure is symmetric with respect to reflection across the ei-axis: B = B\ . 
and F A = A, but -IX + A. 
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Figure 12. Support of the Levy measure for Remark 15.51 showing the 
effect of the exponent B2 on the symmetry group. Here S(jj) = C3. 
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Figure 13. Support of the Levy measure for Remark 15.51 showing the 
effect of the exponent B3 on the symmetry group. Here S(fi) = C\. 

6. Operator self-similar processes 

In this section, we discuss more general operator self-similar processes, whose increments 
need not be independent or stationary. From now on, assume that the operator self-similar 
process X is proper (i.e., for every t > the smallest hyperplane supporting the distribution 
of X(t) equals M. d ), stochastically continuous, and X(0) = 0. Under these assumptions, 
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the real parts of eigenvalues of the exponent B are positive [13j Theorem 4]. Denote by 
<S(X) the set of linear operators A in GL(R d ) such that 

(6.1) {AX(t)} t > f = {X(t)} t > . 

The symmetries of X form a compact subgroup of GL(M d ) as long as X is proper. The 
symmetry group can be seen as minimal information about a multidimensional stochastic 
process. Hudson and Mason [13j Theorem 2] proved that 

(6.2) £(X) = B + TS{X), 

where B G £(X) is arbitrary and T«S(X) is the tangent space of <S(X) at the identity. 
Maejima [23] showed that one can always find a commuting exponent B c G £ (X) such 
that AB C = B C A for all A G «S(X). 

A shift is included in the symmetry group S(fJ,) defined in (|2.3|l for operator stable Levy 
processes, since the definition (12. lh also includes a shift. For operator self-similar processes, 
the definition (jl.lj) does not include a shift, so it is natural that the definition (|6.ip for 
the symmetry group <S(X) of an operator self-similar process does not allow a shift. The 
following lemma connects <S(X) with S(fi) in the operator stable case. 

Lemma 6.1. Let X = {X(t)}t>o be a strictly operator stable Levy process with exponent 
B. Suppose that 1 is not an eigenvalue of B. Then <S(X) = <S(/i), where fi = C(X(1)). 

Proof. Since {X(t)} = {AX(t)} if and only if X(l) = AX(1), we have 5(X) = 5 (//), so 
it suffices to show that S(fi) = So(fi) (see definitions (j2.3jl and (12.71) ). Let A G <S(/x), so 
that AX(1) and X(l) — b are identically distributed for some b G M. d . Since the real parts 
of eigenvalues of all exponents of [i are the same (see (27J Corollary 7.2.12]), we may take 
B as a commuting exponent. Then, for every t > we have 

AX (t) = X(t) -tb = t B X(l) -tb = t B (AX(l) + b)-tb 

= At B X{\) + t B b -tb = AX(t) + t B b - tb. 

Thus (t B - t)b = for all t > 0, and since 1 is not an eigenvalue of B, b = 0. Hence 
A G <So(/i). The converse inclusion, D Sq{h), is obvious. □ 

Remark 6.2. Full dimensional operator stable Levy processes, and proper operator self- 
similar processes, form two distinct classes. Neither class is contained in the other. Take 
Z(t) a spherically symmetric Levy process on M. d whose marginals are Cauchy. Then 
b + Z(t) is an operator stable Levy process, but it is not operator self-similar. The process 
X(t) = Z(t p ) for p > 1 is operator self-similar but not Levy. Remark [6 . 6 1 provides examples 
of operator self-similar processes for which none of the one-dimensional distributions are 
operator stable. The process X(t) = vt + Z(t) is a strictly operator stable Levy process 
and also a proper operator self-similar process. If we take fx = C(X(1)) then <S(/i) = 
but <So(m) = consists of the orthogonal transformations that fix the vector v. 

Theorem 15.11 and Remark 15.41 showed how to construct an operator stable Levy process 
with any admissible symmetry group. The group 0^" (and groups conjugated to it) were 
excluded, since they are not maximal (see Section [T]). This raises a question, is it possible 
to have <S(X) = for some operator self-similar (not Levy) processes? The answer is 
affirmative, as shown in the following example. 
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Example 6.3. Consider a complex valued process 

X(t) = t /3 exp(i(9 + logt)) , t>0, 
where (3 > 0, G is a uniform random variable on [0, 2n] and -X"(0) 

{e^X(t)} t > f ={X(t)} t > , 

X as a process in M 2 , 



0. Since for any 4> £ 



cos(0 + logi) 
sin(B + logi) 



X(t) = t 



is a self-similar with index (5 and C <S(X). By (|6.2 
with b = (3 and arbitrary c). If A € <S(X) then 



I and i?2 are exponents of X (B 2 



which implies A G 02 • Thus 

Consider the process {FoX(t)} t >o, 

F X(t) 



AX(1) = X(1) 



0+ C «S(X) C 2 . 



where Fq is the reflexion with respect to the x-axis, 



p cos(9 + \ogt) 
— sin(G + logt) 

If F G <S(X), then for t\ = \ and £2 = e 7 ^ 2 we would have 

(^(l),^!^ 2 )) = (X(l),X(e^ 2 )), 



or 



cos 
— sin 



sin 
cos 



cos 
sin 



— sin 
cos 



This equality written in M 4 means 



(cos 0, — sin 0, — sin 0, — cos 0) = (cos 0, sin 0, — sin 0, cos 0), 

which is impossible since the sum of the first and the fourth random variables on the left 
hand side is 0, while on the right hand side is 2cos0. Hence Fq ^ <S(X), which yields 
S(X)=0+. 

Remark 6.4. Example l6.3l is consistent with the result that symmetry groups of probability 
measures must be maximal |25l Theorem 2], even though 0^" is not a maximal subgroup of 
GL(M 2 ). This is because, for A in <S(X), we not only require AX(t) identically distributed 
with X(t) for a single t > 0, but also that (AX(t±), . . . ,AX(t p )) is identically distributed 
with (X(t\), . . . ,X(t p )) for all finite- dimensional distributions. We say that 0^" acts diag- 
onally in this case, and we identify A with corresponding element of GL(R 2p ) defined by 
(xi, . . . , x p ) 1— ► (Ax 1, . . . , Axp) for x±, . . . , x p G M 2 . In Example 16.31 the diagonal action of 
O2 is a maximal subgroup of GL(M 4 ), see the proof of Theorem 1 in [2oJ. 



The exponents of a proper operator self-similar process are related to the symmetry 
group by (|6.2p . there always exists a commuting exponent, and the eigenvalues of any 
exponent all have positive real part. These were the crucial facts used in the proof of 
Theorem 15.11 Hence we can also characterize the symmetry group of a proper operator 
self-similar process in M 2 in terms of the exponent B in Jordan form. The proof is identical 
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to Theorem 15,11 except that here we cannot exclude the case where <S(/i) is conjugate to 
C?2 ) as explained in Remark 16.41 

Corollary 6.5. Let X = {X(t)}t>o be a proper operator self-similar process in M? with 
an exponent B given in the Jordan form (15.11) . Then the statements (i)-(iii) of Theorem 
\5.1\ hold verbatim after replacing S{n) by <S(X) and including as a possible symmetry 
group in (ii) 

Remark 6.6. As a simple extension of the construction in Remark 15.41 we can obtain an 
operator self-similar process in M 2 with any exponent, and any admissible symmetry group. 
Take X(t) as in Remark 15.41 and let Y(t) = X(T(t)) where T(t) is a self-similar process 
(time change) with T{at) = a p T(t) (e.g., take T(t) = t p ). Then Y(t) is operator self- 
similar with exponent D = pB. This, together with Example 16.31 also shows that <S(X) 
can take every possible form listed in Corollary 16.51 which therefore provides a complete 
characterization in M 2 of the possible symmetries of an o.s.s. process. An interesting 
and useful example of a self-similar process T(t) with Hurst index < < 1, which is 
not infinitely divisible or even Markovian, is given by the first passage or hitting time 
T(t) = inf{u > : D(u) > t} of a stable subordinator D(t) with E(e~ sD ®) = exp^-cts 13 ). 
The process Y(t) = X(T(t)) has densities h(x,t) that solve the space-time fractional 
multiscaling diffusion equation 

^ ' =Lh(x,t) 

where L is the generator of the operator stable semigroup, see for example pHl ESI E5] ■ This 
fractional diffusion equation models contaminant transport in heterogeneous porous media, 
and the process Y(t) represents the path of a randomly selected contaminant particle. The 
order of the time fractional derivative controls particle retention (sticking or trapping) 
while the exponent of the operator stable process codes the anomalous super diffusion 
caused by long particle jumps. Also, the inverse process T(t) is constant on intervals 
corresponding to jumps of the stable subordinator D(t), the length of which is determined 
by the stable index (3. Note that the time change need not be independent of the outer 
process [3j[33]. Methods for simulating these non-Markovian subordinated processes have 
recently been developed by Magdziarz and Weron [20] and Zhang et al. |46| . 

Remark 6.7. Any proper operator self-similar process can be transformed to one in which 
the exponent takes the Jordan form (|5.ip by a simple change of basis. Corollary 16.51 
shows that the Jordan basis renders the symmetries orthogonal, and then the Jordan form 
determines which symmetries are possible. 

References 

[1] S0ren Asmussen and Jan Rosiriski. Approximations of small jumps of Levy processes with a view 
towards simulation. J. Appl. Probab., 38(2):482-493, 2001. 

[2] William Barker and Roger Howe. Continuous Symmetry: From Euclid to Klein. American Mathemat- 
ical Society, 2007. 

[3] P. Becker-Kern, M.M. Meerschaert and H.P. Schemer (2004) Limit theorems for coupled continuous 

time random walks. The Annals of Probability 32, No. IB, 730-756. 
[4] D. Benson, S. Wheatcraft and M. Meerschaert (2000) Application of a fractional advection-dispersion 

equation. Water Resour. Res. 36, 1403-1412. 
[5] D. Benson, R. Schumer, M. Meerschaert and S. Wheatcraft (2001) Fractional dispersion, Levy motions, 

and the MADE tracer tests. Transport in Porous Media 42, 211-240. 



28 SERGE COHEN, MARK M. MEERSCHAERT, AND JAN ROSINSKI 

[6] P. Billingsley (1966) Convergence of types in k-space. Z. Wahrsch. Verw. Geb. 5, 175-179. 

[7] R. M. Blumenthal and R. K. Getoor, The dimension of the set of zeros and the graph of a symmetric 

stable process. Illinois J. Math. 6 1962 308-316. 
[8] S. Cohen, C. Lacaux, and M. Ledoux. (2008) A general framework for simulation of fractional fields. 

Stochastic Process. Appl. 118(9), 1489-1517. 
[9] S. Cohen and Rosihski. Gaussian approximation of multivariate Levy processes with applications to 

simulation of tempered stable processes. Bernoulli, 13(1):195-210, 2007. 
[10] Rama Cont and Peter Tankov. Financial modelling with jump processes. Chapman & Hall/CRC, Boca 

Raton, Florida, 2004. 

[11] P. Embrechts and M. Maejima (2002) Self-similar Processes. Princeton University Press. 

[12] Franklin A. Graybill. Matrices with applications in statistics. Wadsworth Statistics/Probability Series. 
Wadsworth Advanced Books and Software, Belmont, Calif., second edition, 1983. 

[13] William N. Hudson and J. David Mason. Operator-self-similar processes in a finite-dimensional space. 
Trans. Amer. Math. Soc, 273(l):281-297, 1982. 

[14] H.E. Hurst, R.P. Black, and Y.M. Simaika (1965) Long-term Storage: An Experimental Study, Con- 
stable, London. 

[15] Aleksander Janicki and Aleksander Weron. Simulation and Chaotic Behavior of a- Stable Stochastic 
Processes. Monographs and Textbooks in Pure and Applied Mathematics. Marcel Dekker Inc. New 
York, 1994. 

[16] Zbigniew J. Jurek and J. David Mason. Operator-Limit Distributions in Probability Theory. Wiley 
Series in Probability and Mathematical Statistics. John Wiley & Sons Inc., New York, 1993. 

[17] Olav Kallenberg. Foundations of Modem Probability. Probability and its Applications (New York). 
Springer- Verlag, New York, second edition, 2002. 

[18] Peter E. Kloeden and Eckhard Platen. Numerical solution of stochastic differential equations. Appli- 
cations of Mathematics. Springer- Verlag, New York, 1992. 

[19] Celine Lacaux. Series representation and simulation of multifractional Levy motions. Adv. in Appl. 
Probab., 36(1):171-197, 2004. 

[20] Magdziarz, M., and A. Weron, Competition between subdiffusion and Levy flights: A Monte Carlo 
approach, Phys. Rev. E, 75, 056702, 2007. 

[21] Maejima, M. and J. D. Mason (1994) Operator-self-similar stable processes. Stoch. Proc. Appl. 54, 
139-163. 

[22] Makoto Maejima. Operator-stable processes and operator fractional stable motions. Probab. Math. 

Statist., 15:449-460, 1995. Dedicated to the memory of Jerzy Neyman. 
[23] Makoto Maejima. Norming operators for operator-self-similar processes. Trends Math., Birkhauser, 

Boston 287-295, 1998. 

[24] M.M. Meerschaert and H.P. Scheffler, One dimensional marginals of operator stable laws and their 

domains of attraction, Publ. Math. Debrecen, 55(3-4), 487-499, 1999. 
[25] Mark M. Meerschaert and J. A. Veeh. Symmetry groups in d-space. Statistics & Probability Letters 

22:1-6, 1995. 

[26] Mark M. Meerschaert and Hans-Peter Scheffler. Sample cross-correlations for moving averages with 
regularly varying tails. Journal of Time Series Analysis, 22(4):481-492, 2001. 

[27] Mark M. Meerschaert and Hans-Peter Scheffler. Limit Distributions for Sums of Independent Random 
Vectors: Heavy Tails in Theory and Practice. Wiley Series in Probability and Statistics. John Wiley 
& Sons Inc., New York, 2001. 

[28] Meerschaert, M., D. Benson and B. Baeumer (2001) Operator Levy motion and multiscaling anoma- 
lous diffusion. Phys. Rev. E 63, 1112-1117. 

[29] M.M. Meerschaert, D.A. Benson, H.P. Scheffler and B. Baeumer (2002) Stochastic solution of space- 
time fractional diffusion equations. Phys. Rev. E 65, 1103-1106. 

[30] M.M. Meerschaert and H.P. Scheffler, Portfolio modeling with heavy tailed random vectors. Handbook 
of Heavy-Tailed Distributions in Finance, 595-640, S. T. Rachev, Ed., Elsevier North-Holland, New 
York, 2003. 

[31] Mark M. Meerschaert and Yimin Xiao. Dimension results for sample paths of operator stable Levy 
processes. Stochastic Processes and Their Applications, 115(l):55-75, 2005. 

[32] M.M. Meerschaert, E. Scalas, Coupled continuous time random walks in finance. Physica A: Statistical 
Mechanics and Its Applications, 370, 114-118, 2006. 



MODELING AND SIMULATION WITH OPERATOR SCALING 



29 



[33] Meerschaert, M.M. and H.-P. Scheffler (2008) Triangular array limits for continuous time random 

walks. Stock. Proc. Appl., 118, 1606-1633. 
[34] W.E. Pruitt and S.J. Taylor (1969) Sample path properties of processes with stable components. Z. 

Wahrsch. verw. Geb. 12, 267-289. 
[35] Rachev, S. and S. Mittnik (2000) Stable Paretian Models in Finance, Wiley, Chichester. 
[36] B.S. Rajput and J. Rosihski. Spectral representations of infinitely divisible processes. Probab. Th. Rel. 

Fields, 82: 451-487, 1989. 

[37] D.M. Reeves, D.A. Benson, M.M. Meerschaert, H.P. Scheffler, Transport of Conservative Solutes in 
Simulated Fracture Networks 2. Ensemble Solute Transport and the Correspondence to Operator- 
Stable Limit Distributions, Water Resources Research, 44 (2008), W05410. 

[38] J. Rosihski. On series representations of infinitely divisible random vectors. Ann. Probab., 82: 405-430, 
1990. 

[39] Jan Rosihski. Series representations of Levy processes from the perspective of point processes. In Levy 

processes, pages 401-415. Birkhauser Boston, Boston, MA, 2001. 
[40] Jan Rosihski. Tempering stable processes. Bernoulli, 13:195-210, 2007. 

[41] Ken-Iti Sato. Strictly operator-stable distributions. Journal of Multivariate Analysis, 22:278-285, 
1987. 

[42] Ken-iti Sato. Levy Processes and Infinitely Divisible Distributions, volume 68 of Cambridge Studies 
in Advanced Mathematics. Cambridge University Press, Cambridge, 1999. Translated from the 1990 
Japanese original, Revised by the author. 

[43] E. Scalas, R. Gorenflo, F. Mainardi, and M.M. Meerschaert, Speculative option valuation and the 
fractional diffusion equation, Fractional Derivatives and Their Applications, 265-274, A. Le Mehaute, 
J. A. Tenreiro Machado, J. C. Trigeassou and J. Sabatier, Eds. (2005), Ubooks, Germany. 

[44] Oleg Sheluhin, Sergey Smolskiy, Andrew Osin, Self-Similar Processes in Telecommunications, Wiley, 
New York, 2007. 

[45] Y. Zhang, D.A. Benson, M.M. Meerschaert, E. M. LaBolle, and H.P. Scheffler. Random walk approx- 
imation of fractional-order multiscaling anomalous diffusion. Physical Review E, 74(2):6706-6715, 
2006. 

[46] Y. Zhang, M.M. Meerschaert, B. Baeumer (2008) Particle tracking for time-fractional diffusion, Phys- 
ical Review E, 78(3), 036705. 

Serge Cohen, Universite de Toulouse, Universite Paul Sabatier, Institut de Mathema- 
tiques de Toulouse. F-31062 Toulouse. France. 
E-mail address: Serge.Cohen@math.univ-toulouse.fr 
URL: http: //www. math. univ-toulouse . f r/~cohen/ 

Mark M. Meerschaert, Department of Mathematics, University of Nevada, Reno NV 
89557 USA 

E-mail address: mcubed@stt.msu.edu 
URL: http : // www . stt . msu . edu/ ~mcubed 

Jan Rosinski, Department of Mathematics, University of Tennessee, Knoxville, TN 
37996, USA. 

E-mail address: rosinski@math.utk.edu 
URL: http: //www. math. utk. edu/~rosinski/ 



